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Abstract 

We present first results for B c spectroscopy using Lattice Non-Relativistic 
QCD (NRQCD). For the NRQCD action the leading order spin-dependent 
and next to leading order spin-independent interactions have been included 
with tadpole-improved coefficients. We use multi-exponential fits to mul- 
tiple correlation functions to extract ground and excited S states and give 
accurate values for the S state hyperfine splitting and the P state (B**) 
fine structure, including the effects of 1 P\/^P\ mixing. 
PACS : 12.38.Gc, 12.39.Hg, 14.40.Lb, 14.40.Nd 



1 Introduction 



There is much current interest in the spectrum of mixed bound states of bottom 
and charm quarks with the recent appearance of experimental candidates Lat- 
tice QCD is the best method of calculating this spectrum from first principles, 
although attention must be paid to the control of systematic errors from a variety 
of sources. We have recently given accurate results for bottomonium and char- 
monium spectroscopy on the lattice 0, [!|, exploiting the non-relativistic nature 
of these systems by using a non-relativistic effective theory, NRQCD |[|. This can 
be systematically matched to full QCD, order by order in a s and v 2 /c 2 , where v 
is the velocity of the heavy quark in the bound state. Here we present results for 
B c spectroscopy using the same techniques and use our previous results to give 
good indications of the size and direction of systematic errors. We compare our 
results to recent calculations in potential models 0. 
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2 Lattice NRQCD 



The starting point of NRQCD is to expand the original QCD Lagrangian in 
powers of v 2 , the typical quark velocity in a bound state. For the J/ty system 
v 2 ~ 0.3 and for T v 2 ~ 0.1. For B c v 2 ~ 0.15 for the single particle with mass 
equal to the reduced mass of the b, c system, but the kinetic energy will be shared 
unequally between the b and c. v 2 ~ 0.5 for the c quark so relativistic corrections 
will be more important than for cc systems. NRQCD enables these corrections 
to be added systematically. 

The action used here is the same one used in refs. and 0. Relativistic 
corrections 0(Mqv a ) have been included for both quarks - errors are then at 
the level of 0(M c v®). Other sources of systematic error are discretisation errors 
and errors from the absence of virtual quark loops because we use quenched 
configurations generated with the standard plaquette action. Finite volume errors 
should be negligible because of the relatively small size of heavy-heavy systems. 

To calculate masses for be bound states we define b and c quark Green func- 
tions on the lattice. The NRQCD Lagrangian involves a simple difference equa- 
tion in the temporal direction, which allows the evolution of the quark Green 
function as an initial value problem. We start on the first time slice with 

and then continue to evolve in the temporal direction using 

G t+1 = (l-p) n ut(l-^) n (l-a5H)G t (t > 1). (2) 



2n J * V In 

Hq is the lowest order piece of the non-relativistic Hamiltonian, i.e. the kinetic 
energy operator. On the lattice, this is 

A( 2 ) 

Ho = "2Mf (3) 
The higher order terms in the Hamiltonian that we have included are 

SH = -Cx^^ + Cs-^fA-E-E- A) 
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The first two terms in 5H are spin-independent relativistic corrections and the 
next two are spin- dependent terms which contribute the leading order pieces to 
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the P and S spin splittings respectively. The last two terms come from discretisa- 
tion corrections to the lattice Laplacian and the lattice time derivative. A is the 
symmetric lattice derivative, A^ 2 -* is the lattice form of the Laplacian and A*- 4 - 1 is 
a lattice version of the continuum operator J2Df. We use the standard traceless 
cloverleaf operators for the chromo-electric and magnetic fields, E and B. The 
parameter n is introduced to remove instabilities in the heavy quark propagator 
caused by the highest momentum modes of the theory [Q. We require n > 3 /Ma. 
To calculate we must fix both the bare coupling and the bare quark masses. We 
used 200 quenched gluon field configurations of size 12 3 x 24 at (3 = 6/g 2 = 
5.7 generously supplied by the UKQCD collaboration |J, and fixed to Coulomb 
gauge, using a Fourier accelerated steepest descents algorithm 0. The bare 
masses of the c and b quarks were chosen from fits to the kinetic mass of the r] c 
and the T 0, |. At (3 = 5.7 we used a mass for the c quark in lattice units of 
0.8 with n = 4, and a mass for the b quark of 3.15 with n = 2. 

The coupling constants c* appearing in equation can be calculated by 
matching NRQCD to full QCD [|| [IU]. At tree level all the coefficients are one. 



The largest radiative corrections are believed to be tadpole contributions which 



can be removed by redefining the gluon fields U : 11 



U,(x) - ^ (5) 

with uq the fourth root of the plaquette (at (3=5.7 we use u Q = 0.861). Since 
the cloverleaf expression involves the evaluation of a plaquette this renormaliza- 
tion will have a large effect on E and B fields and thereby on spin-dependent 
splittings With the dominant tadpole contributions thus removed, we use the 
tree level values for the q's. 

Given the quark propagators in equation (|2]) it is relatively straightforward to 
combine them appropriately to form meson propagators with specific quantum 
numbers. This is outlined in ||. Here we must combine a b quark propagator 
with a c propagator, but the methods are identical. We use 'smearing functions' 
as sources for the b quark propagator and a delta function source for the c quark. 
Various different smearing functions are included to study p states and excited s 
states. It is more efficient, from equation (0), to have numerous b quark propa- 
gators with a small value of n and only one c quark propagator. It also enabled 
us to simultaneously calculate the bb spectrum on these lattices, and this will be 
reported separately |§ . The radius of the smearing functions (taken as wavefunc- 
tions from a 1/r potential) was kept fixed at an optimal value for bb correlation 
functions. This meant that the effective mass in the B c correlation functions did 
not plateau quite as early as it would have done with more optimal smearing but 
it did not significantly detract from the accuracy of our results. 

Local meson operators are tabulated in M. Using the notation 2S+1 Lj, we 
have looked at meson propagators for the following states: 1 S , 3 Si, 1 Pi, 3 P , 3 Pi 
and 3 P2 for both the E and T representation. Since C is not a good quantum 
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number for the B c system the X P\ and 3 Pi mesons will mix and so in addition 
we calculated the cross-correlation between these two. For the s states, smearing 
functions both for the ground and first radially excited state were used as well as 
a local S function. From this all possible combinations of smearing at the source 
with a local sink were formed. For the p states only the ground state smearing 
function was used at the source. We calculated the dispersion relation for the 1 S 
(B c ) by looking at the meson propagator for small non-zero momenta. Because 
of the relatively small size of these systems it is possible to use more than one 
starting site for the mesons. We used 8 different spatial origins and 2 different 
starting times to increase our statistics, but we bin correlation functions over 
each configuration. 

3 Simulation results 

For bb and cc the correlation functions are explicitly real because of charge con- 
jugation symmetry - for be they are not. We find, however, that the imaginary 
parts of all our correlation functions are very small and much smaller than the 
real parts. We fit simultaneous multi-exponential fits to the real parts of several 
correlation functions with different sources and local sinks. This is the 'row' fit 
ofrefs. g§ : 

N exp 

G meson (n sc ,loc;t) = b(n sc ,k)e~ Ek4 . (6) 

fc=l 

n sc denotes the source smearing, loc denotes a local sink. It enables us to ex- 
tract masses in lattice units for ground and excited s states. The p states are 
much noisier and for them we have used single exponential fits to the correlation 
function with ground state smearing at both source and sink. For spin splittings 
it is most accurate to perform a single exponential fit to a ratio of correlation 
functions. We generate a bootstrap ensemble of ratios (including both real and 
imaginary parts of the correlation functions) and then fit to the real part. 

To find masses for the eigenvectors of the l P\ and 3 Pi mixing matrix we fit a 
2-exponential form to the 2x2 matrix formed from these correlation functions 
smeared at source and sink. Figure 1 shows effective amplitude plots and fits 
for the four correlation functions in the 2x2 matrix. The physical states are 
called the 1 + and the 1 + . To calculate the splitting between these two states 
and the mixing angle with the l Pi/ 3 Pi basis accurately, we did matrix fits to an 
ensemble of matrices generated by bootstrap. This allows the correlation between 
the states on a given configuration to be taken into account. We checked these fits 
against those obtained by diagonalising the matrix time-slice by time-slice and 
fitting single exponential forms to the eigenvalues. There was good agreement 
between the two methods. We did not find it possible to extract the masses of the 
1 + and 1 + purely from the 1 P\ and 3 Pi correlation functions (i.e. the diagonal 
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Figure 1: Effective amplitudes (obtained by multiplying the correlation function 
by exp(M groun( }t)) for the four correlation functions in the 2x2 1 P 1 / 3 P 1 mixing 
matrix. 



terms of the matrix) in the absence of the cross-correlation terms. Although 
theoretically possible this obviously requires very high statistics or a more finely 
grained lattice in time. 

Table 1 shows our results for masses and splittings in lattice units. The errors 
quoted are statistical errors only. 

To convert the masses and splittings in lattice units to physical results, we 
need a value for a -1 , the inverse lattice spacing. Heavy-heavy systems provide a 
particularly good way of fixing a -1 because the spin-averaged IP-IS splitting is 
insensitive to the other parameter in the theory, i.e. the heavy quark mass, and 
also to any systematic errors in spin-dependent terms. Thus a~ l can be fixed from 
a comparison of either the cc or the bb spectrum to experiment. Experimentally 
the spin averaged splitting is 452 MeV from the bb spectrum and almost the same, 
458 MeV, from the cc spectrum. Here the spin-averaged IP mass is taken as the 
average of the Xb states, since only a preliminary value for the mass of the h c is 
known and no hb has been seen [[L2| . Potential models also lead to the belief that 
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Simulation Results 


l^So 


0.6052(8) 


1 3 S\ 


0.6353(9) 


2 1 S'o 


1.12(6) 


2 3 5i 


1.14(6) 


1 3 P 2 


0.986(7) 


1 3 P 


0.944(7) 


1 1+ 


0.956(6) 


1 1+' 


0.973(7) 


3C 1 Q 


0.0307(2) 


3 P 2 - 3 Po 


0.045(4) 


3 P 2 - 1 + 


0.023(2) 


1+' - 1 + 


0.0172(25) 


tan(9) 


0.66(3) 


m Bc 


0.2118(4) 



Table 1: Fitted dimensionless energies for be states and mixing angle for the 1 + 
state. 9 is defined so that |1 + ) = sin(9)\ 1 Pi) + cos(9)\ 3 Pi). The final entry is the 
wavefunction at the origin for the B c . 



the l P\ state should lie approximately at the spin-average of the 3 P states, and 
so does not need to be included. 

Unfortunately, in the quenched approximation, the value obtained for a -1 
depends on the experimental quantity you use to fix it, and in particular, the 
momentum scales important to that quantity. This is because the effective cou- 
pling constant does not run correctly between different momentum scales. In 
general you expect that a quantity which is more sensitive to large momenta will 
give a higher value for a -1 . This is what we find on comparing the bb and cc 
spectrum on quenched lattices. At (3 = 5.7 = 1.20(4) GeV [|3| and a^- 1 = 
1.35(4) GeV [f|. These results use the spin-average of all IP states [13| because 
we find that the l P\ does not lie at the spin average of the 3 P states. This is 
likely to be, at least partly, a discretisation error ||. 

This discrepancy in values obtained for a~ l from the quenched approximation 
means that the only way to obtain results in good agreement with experiment 
for a particular set of hadrons is to fix a~ x separately for that set using one 
experimental result for a 'typical hadron' from the set. For be systems there are 
no experimental results so we must fix a^ 1 on the basis of the numbers for bb and 
cc. We assume that the spin-averaged IP-IS splitting for be will have a value 
between that for bb and that for cc (which hardly differ, as above). We must now 
define the spin-averaged P state from our results to include all IP states since 
there is no reason for either of the 1 + or the 1 + to lie at the spin-average of the 
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other states. This gives a -1 = 1.32(4) GeV, as expected between the values from 
cc and bb. The error quoted is purely statistical. 

The difference in a -1 between these different systems leads to the effect that 
the quark mass in physical units will differ between be and either bb or cc. This 



problem has already been pointed out for heavy-light systems [II] and means that 
again, in principle, one should fix the quark mass, using experiment, separately 
within each system in the quenched approximation. For the B c system, this is 
not possible, so final meson masses will have an additional uncertainty associated 
with this a -1 effect. It will be important particularly for spin splittings since these 
are strongly dependent on the quark mass. We have fixed the bare quark masses 
in lattice units from results for bb and cc. Another reasonable possibility is that 
we should readjust the bare quark masses in be so that they agree in physical 
units with those of bb and cc. We have estimated what changes this might give 
where possible. 

Using the value for a~ l above, we can convert meson masses to physical units. 
We still need to set a zero of energy since we have removed it from our Hamil- 
tonian. We do this by calculating B c correlation functions at small but non-zero 
values of the meson momentum, P, and fitting to 

E E P2 (7) 

M^in is the kinetic (absolute) mass of the B c in lattice units. We find a kinetic 
mass of 4.76(2), which, using a -1 as above, gives a mass in physical units of 
6.28(20) GeV, where the error is dominated by the uncertainty in a -1 . Errors in 
the b and c quark masses discussed above actually cancel in the kinetic mass of 
the B c , since to match the bare physical masses from bb and cc we in fact have 
to adjust the c mass down and the b mass up by the same amount. 

The difference between Eq and Mki n represents the shift of the zero of energy. 
This can be calculated per quark in perturbation theory and agreement with 
lattice results is very good for the T spectrum 0. For c quarks the perturbative 
results are not very reliable because the important scales for a s are very low. 
However, we can still check that the non-perturbative energy shift for the B c is 
the average of those for T and ^. This should be the case if no non-perturbative 
effects specific to the mesons appear, and the results in Table 2 show that this 
is true. It is a requirement for the non-relativistic theory to make sense and is 
apparently not obeyed for heavy quark actions based on the Wilson action in 
current use[^ . 

Table 3 gives results for the meson splittings with the ground state (B c ) in 
GeV and the spectrum is plotted in Figure 2 (using the B c mass of 6.28 GeV), 
with a comparison to results from a recent potential model analysis fl5|. 
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E 


Mkin 


shift 


T 


0.5030(5) 


6.97(8) 


6.47(8) 


Vc 


0.618(1) 


2.430(6) 


1.812(6) 


B c 


0.6052(8) 


4.76(2) 


4.15(2) 



Table 2: Results for static and kinetic masses and the difference between them 
for heavy-heavy mesons, using NRQCD at (3 =5.7. Values are in lattice units. 





Simulation Results [GeV] 




-B c 


0.0405(3) 


2% 


-B c 


0.68(8) 


2 3 S 1 


- B c 


0.71(8) 


l 3 Po 


-B c 


0.447(9) 


1 1+ 


-B c 


0.463(8) 


1 1+' 


-B c 


0.485(8) 


1 3 P 2 


-B c 


0.503(8) 



Table 3: NRQCD results for splittings of bc states with B c for a 1 = 1.32 GeV. 
Errors shown are statistical only. 

4 Discussion 

Figure 2 shows the lattice spectrum compared to recent potential model results 
0. The lattice results have systematic errors resulting from higher order rela- 
tivistic corrections to the action, discretisation effects and errors from using the 
quenched approximation. We will attempt to quantify them below. The potential 
model results have systematic errors also from various sources including varia- 
tions in the potential itself (compare, for example flB]] ). The lattice results have 
the advantage that the errors there can, and will, be systematically removed. 

The main source of error for our lattice calculations is in higher order rel- 
ativistic corrections to the charm quark propagator. The missing terms are at 
order Mt> 6 and include D /16M|, relevant for spin-independent splittings. A 
perturbative estimate of the size of this correction would be an energy shift of 
< p 6 > /16M^. A naive evaluation of this expectation value in a lattice potential 
model for B c gives 100 MeV. However, there are several terms at this order that 
should be included and we would expect some cancellation between them. A 
similar analysis of terms of order Mv 4 , which are included here, shows an almost 



complete cancellation between the D 4 correction and the Darwin term |T7| . Our 
comparison of lattice results with and without these Mv 4 terms shows agreement 
with this analysis @, giving us confidence that we can estimate these correc- 
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tions. The conclusion is then that terms of order Mv 6 not included could shift 
our masses by 100 MeV, but by cancellation it could be less than this. 

The leading discretisation error is an 0(a 2 ) error from using gluon field config- 
urations with the simple plaquette action. Again we can estimate perturbatively 



the shift induced by this correction |I8|. The shift is proportional to the square 
of the wave function at the origin, giving 15 MeV for the B c , and zero for states 
with no wave function at the origin such as p states. The quark propagators are 
already correct through 0(a 2 ) and residual a 4 errors should be smaller than this. 

The errors from the quenched approximation should be between those for bb 
and cc and can be estimated from the differences we find in simulating those 
systems with and without dynamical fermions. One source of error, that of 
uncertainties in the quark masses, has been discussed above. 

From Figure 2 the absolute mass of the B c which we obtain agrees well with 
potential model estimates. Since it is close to the average of the r\\, and r\ c masses 
this is not surprising. The spin-averaged IP-IS splitting also agrees by design, 
since we fixed it to be the same as that for bb and cc and this would come out 
of a potential model too. The possible 100 MeV (i.e. 25%) systematic error in 
the IP-IS splitting from higher order relativistic corrections should be borne in 
mind. The corresponding systematic error for bb and cc is 6 MeV and 40 MeV 
respectively. 

The other spin-independent splitting is that between the 2S and IS states. All 
our quenched calculations || |], § give a ratio for the 2S-1S splitting over IP-IS 
splitting which is too large. This is expected because the IS state suffers a bigger 
shift downwards under quenching than the other states. We would therefore 
expect our 2S states to appear too high. For bb we find the 2S-1S/1P-1S ratio 
(using the 3 S 1 for S) to be 1.40(3) at (3 = 6.0 |L9| and 1.4(2) at (3 = 5.7 [§ 
compared to the experimental value of 1.28 and for cc we get 1.4(2) |§ against 
the experimental value of 1.38. So, we would expect that unquenching would 
correct this ratio for B c by roughly 10%. Since we use the IP-IS splitting to fix 
a -1 , this means that our 2S-1S splitting would be reduced by 10%, i.e. our 2S 
level would fall by 70MeV. This brings it rather closer to the potential model 
prediction, which presumably already incorporates some effects from dynamical 
fermions through the phenomenological form of the potential. Large relativistic 
corrections for B c could distort this picture somewhat. 

The spin splittings are perhaps of more immediate interest and particularly 
that between the B* and the B c , since these states will be the first to be studied 
experimentally in detail. We find 40 MeV for this splitting in this simulation. 
Again we expect a sizeable error from the quenched approximation. For bb we 
have results both for dynamical flavors, Nf = and 2 and find that the hyperfine 
splitting changes by a factor of 40 % on unquenching. For cc the hyperfine 
splitting is known and our quenched simulation gives a result 20 % too low. So 
we conclude that the splitting between the B* and B c might reasonably change 
to 55 MeV on unquenching. It will increase by an additional 10% to 60 MeV 
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if the charm quark mass in the B c is reduced so that it matches the bare value 
in physical units for the J/^f. Both of these effects again tend to improve the 
agreement with the potential model (perturbative) value. 

For the p fine structure quenching effects may not be very large because of the 
small wavefunction at the origin. More serious are discretisation errors since the p 
fine structure is the result of a balance between short and long range effects. From 
work on bb [§, §] and cc || we find a tendency for lattice NRQCD calculations 
to underestimate the overall size of the splittings [M( 3 P 2 ) — M( 3 P )], a larger 
effect for bb at (3 = 5.7 than for cc. We believe that this is because the short- 
range components of the spin-orbit and tensor forces are underestimated at this 
lattice spacing. If so, the 3 P 2 — 3 Po splitting will also be underestimated here, 
from the same effect. The result for the splitting would additionally increase if 
the charm quark mass were reduced. Our value for this splitting is 60(5) MeV, 
already larger than some potential model results. It seems likely to us that the 
experimental value will exceed 60 MeV. 

Another result from our bb and cc work was the ratio of fine structure splittings 
[M(x 2 ) — M(xi)]/[M(xi) — M(xo)]- This exceeded the experimental values and 
was another indicator that the long-range spin-orbit force has undue dominance 
at this lattice spacing. For the perturbative potential model of ref. H this ratio 
comes out below experiment for charmonium (although not for bottomonium) , 
indicating that perhaps the long-range spin-orbit term is not large enough. 

For B c the xi state mixes with the l P\. Our results show such a large amount 
of mixing that the physical states 1 + and 1 + are very close to the jj coupled 
states expected in the M b — > oo limit. For these states we couple L + s c = J c and 
J c + si = J ■ The J c = 1/2 state has a mixing angle 9 with the LS coupled basis 
of tan(0) = l/y(2) = 0.71. In the Mb — ► oo limit the 1 + state becomes the J c 
= 1/2 state, degenerate with the Xo, an d the 1 + ', the J c = 3/2 state, degenerate 
with the Xi- Our result of tan(#) = 0.66 is quite different to the (perturbative) 
potential model || in which the mixing (viewed from the LS coupled basis) is 
very small. 

The true result probably lies somewhere in between. In the LS basis and 
in potential model language, the off-diagonal terms of the mixing matrix are 
provided purely by the long-range spin-orbit potential. The short-range spin- 
orbit and tensor terms contribute only to the diagonal pieces (as does the spin- 
spin term, although we expect it to vanish for p states). For our calculation 
then the off-diagonal terms are too large compared to the on-diagonal, and we 
overestimate the mixing. For the perturbative potential models, the opposite 
may be true. The mixing is clearly very sensitive to these effects. 

We can also estimate the decay constant fs c from our result for the wave- 
function at the origin. Using the standard formula f 2 = 12|?/>(0)| 2 /M, we obtain 
fs c = 440(20) MeV, where the error quoted is dominated by the uncertainty in 
a -1 . This result is only valid to leading order in the inverse quark masses and 
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B c B* Pstates 



Figure 2: NRQCD simulation results for the spectrum of the B c system using 
an inverse lattice spacing of 1.32 GeV. Error bars are shown where visible and 
only indicate statistical uncertainties. Dashed lines show results from a recent 
potential model calculation [|5| 

a s . Relativistic 1/M% corrections should be applied as well as a perturbative 
renormalisation factor before a useful continuum value can be quoted ||20|| . 

5 Conclusions 

This represents a first calculation of the be spectrum on the lattice. We use 
NRQCD and include the leading relativistic and discretisation corrections with 
tadpole- improved coefficients. We give a spectrum including radially excited s 
states as well as p fine structure taking account of mixing between the J = 1 
states. Agreement with potential model results is surprisingly good, given that 
they include, at least explicitly, no relativistic corrections and the velocity of the 
charm quark within a B c is actually higher than for charmonium. Analysis of the 
systematic errors in the lattice calculation tends to improve the agreement with 
potential models except for the overall scale of the p fine structure, M( 3 P 2 ) — 
M( 3 P ). We believe this is underestimated at present. Our result for the mixed 
J = 1 states show strong mixing in the LS coupled basis so that physical states 
are close to the jj coupled limit. We believe that the improvement of systematic 
errors from discretisation will tend to reduce this mixing. Future calculations of 
the spectrum will work at smaller lattice spacings with relativistic charm quarks. 
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